Critical behavior of a colloid-polymer mixture confined between walls 
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We investigate the influence of confinement on phase separation in colloid-polymer mixtures. To 
describe the particle interactions, the colloid-polymer model of Asakura and Oosawa [J. Chem. 
Phys. 22, 1255 (1954)] is used. Grand canonical Monte Carlo simulations are then applied to this 
model confined between two parallel hard walls, separated by a distance D = 5 colloid diameters. 
We focus on the critical regime of the phase separation and look for signs of crossover from three- 
dimensional (3D) Ising to two-dimensional (2D) Ising universality. To extract the critical behavior, 
finite size scaling techniques are used, including the recently proposed algorithm of Kim et al. 
[Phys. Rev. Lett. 91, 065701 (2003)]. Our results point to "effective" critical exponents that differ 
profoundly from 3D Ising values, and that are already very close to 2D Ising values. In particular, 
we observe that the critical exponent /3 of the order parameter in the confined system is smaller 
than in 3D bulk, yielding a "flatter" binodal. Our results also show an increase in the critical 
colloid packing fraction in the confined system with respect to the bulk. The latter seems consistent 
with theoretical expectations, although subtleties due to singularities in the critical behavior of the 
coexistence diameter cannot be ruled out. 

PACS numbers: 05.70.Jk, 64.60.Fr, 64.70.Fx 



I. INTRODUCTION 



The current technological demand for the production 
of nanoscopic devices [J. H H H renewed the inter- 
est in understanding the phase behavior of fluid sys- 
tems confined in pores of nanoscopic linear dimensions 
0,0- In addition, porous materials with pore widths 
less than 50 nm are widely used in the chemical, oil and 
gas, food, and pharmaceutical industries, for applications 
such as mixture separation, pollution control, and as cat- 
alysts nuns. However, many such applications rely 
largely on empirical knowledge, since the theory-based 
understan ding of confined fluids is still rather incomplete 
[1 IH 13 El El 111 Even the basic phenomenon of 
"capillary condensation" of undersaturated gases in cap- 
illaries, described already in the 19 th century 01 1 still 

an- 
and 

computer simulations J^ilillMIMIIlElE El 

Regarding confined binary mixtures, there is a 
close analogy between the phase behavior of confined one- 
component fluids and the preferential adsorption of one 
of the components of the mixture to the walls. The misci- 
bility gap of the mixture corresponds to the coexistence 
curve (or binodal) that describes the phase separation 
between liquid and gas in simple fluids, and numerous 
theoretical and simulational studies have addressed the 
phase behavior of binary mixtures in cylindrical pores or 

slit pores [H m m Hi m m El El El El E1E3- 

Depending on the details of the wall-particle interactions 
in relation to the interactions among the fluid particles 
of a binary (AB) mixture, it is clear that there can be 



forms the subject of longstandin 
alytical theory 0, ILsL [l 



meting investigations b 




either the A-component or the B-component attracted 
to the walls, apart from the very special case of "neutral 
walls" which produce confinement only |4g, 1491 150| . Sim- 
ilarly, in a one-c omp onent liquid-gas system "capillary 
evaporation" |4g, |5l|, IHl^I can occur for repulsive wall- 
particle interactions. While for the liquid-gas transition 
and the demixing transition of binary fluids the "order 
parameter" of the transition is a simple scalar, i.e. the 
transition belongs to the "universality class" of the Ising 
(lattice gas) model [HI] , related phenomena occur in sys- 
tems with more complex ordering, e.g. confined liquid 
crystals where "capillary nematization" may occur [55j. 

Understanding nanoscopic confinement of fluids con- 
sisting of small molecules is difficult because the lateral 
variation of the wall potential, due to wall roughness or 
even atomistic corrugation |56| of the wall potential, may 
have drastic__effects on the phase behavior of the con- 
fined fluid [57|. In this respect, colloidal systems due 
to the mesoscopic size of the colloidal particles pose dis- 
tinct advantages: atomistic corrugation of the confining 
wall potentials can safely be neglected, and there is a 
great freedom in preparing systems with suitable interac- 
tions HI IH |H • A particularly suitable class of systems 
are colloid-polymer mixtures, since both bulk phase be- 
havior and the interfaces separating the colloid-rich and 
polymer-rich phases can be studied experimentally in de- 
tail 0, |6l lo f&p. Furthermore, the Asakura- Oosawa 
(AO) model |65l l6q provides a simple theoretical de- 
scription, which seems to capture all the salient features 
of such phase-separating colloid-polymer mixtures, and 
which is well suited to computer simulation investiga- 
tions [Ililillil- 

Therefore, we use this model again in the present pa- 
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per to address the question: how does the critical behav- 
ior near the demixing critical point change due to the 
confinement in slit pores? While on general theoretical 
grounds one expects [II G3 B III El El 111 that in 
the ultimate vicinity asymptotically close to the critical 
point the critical exponents of the two-dimensional (2D) 
Ising model should apply [70| . this limiting asymptotic 
behavior may be hard to observe, and the need arises 
to consider the crossover from three dimensional (3D) 
to two-dimensional critical behavior in such thin films 
HH El EH ■ This crossover behavior in itself is a difficult 
but interesting problem |49|. Previous work considering 
this issue has been restricted to either simple Ising (lat- 
tice gas) models 0, El EH or strictly symmetric poly- 
mer mixtures confined by neutral walls |5(l | . The compli- 
cations due to the asymmetry between liquid(-like) and 
gas(-like) phases in the bulk have not been considered, 
and also the further asymmetry arising from the pref- 
erential adsorption of one species to the walls has been 
disregarded. Due to this preferential adsorption, one ex- 
pects near the critical point of the bulk the formation of 
wetting layers at the walls El El El El E| 

if the width 

of the slit pore becomes very large. 

We note that previous theoretical work on capillary 
condensation (or evaporation) based on density func- 
tional theory (DFT), or other analytical approximations 
for the equation of state, inevitably implies a parabolic 
shape of the binodal near the critical value r]p cr of the 
polymer reservoir packing fraction n^, irrespective of 
whether one considers the bulk mixture or a confined sys- 
tem. Therefore, despite the fact that such theories are 
very powerful away from the critical point, they cannot 
be used to describe the crossover in the critical behavior 
due to confinement. Also previous Gibbs ensemble sim- 
ulations of the AO model confined to slit pores Elj did 
not address this issue since the Gibbs ensemble cannot 
be used to sample the critical point. In the present work, 
we therefore extend the grand canonical techniques used 
in our previous work on the critical behavior of the AO 
model in the bulk [6?l l68| to study the critical behavior 
in confinement. 

The outline of this paper is as follows. In Section ITU 
we briefly review crossover scaling relations that are ex- 
pected to describe the critical behavior of a confined fluid. 
Next, we introduce the AO model and describe our sim- 
ulation method. In Section IlIII we present our results, 
using a finite size scaling analysis of the critical proper- 
ties for a very thin film of thickness D = 5er c , with o~ c 
the diameter of the hard-sphere colloids. We end with a 
summary and conclusion in the last section. 

II. CROSSOVER SCALING 

To define the problem of study more precisely, we recall 
that the width of the binodal, or order parameter, of a 
colloid-polymer mixture is given by 

A = {ri - rf c )/2, (1) 



with n l c (rjc) the colloid packing fraction in the colloidal 
liquid (vapor) phase. In 3D bulk, the order parameter 
close to the critical point is expected to scale as 

A(OO) = Bt£,, <oc = r?pA?p,cr(°°) - 1, (2) 

where B is a (nonuniversal) critical amplitude, and /3 « 
0.326 the (universal) critical exponent of the 3D Ising 
universality class El EH- The symbol (oo) in the 
above emphasizes that 7y p cr is the critical value of the 
polymer reservoir packing fraction appropriate for an in- 
finitely thick film, i.e. a bulk 3D system. 

In a confined thin film of thickness D, however, the 
corresponding relation reads 

A(D) = B(D)t%, t D =rf p /r,l CI {D)-l. (3) 

The critical polymer reservoir packing fraction is thus 
shifted from its bulk value ?7p C1 .(oo) to a new value 

rf pcr (D). In addition, the critical amplitude B(D) de- 
pends on the film thickness D, and the critical exponent 
takes the value of the 2D Ising universality class /?2 = 1/8 
[3 EH- However, as D gets large, the validity of Eq.|j3J) 
is expected to be observable only in an extremely nar- 
row region around tjj = 0. This is recognized when one 
formulates the app ropriate crossover scaling description 

BEEIHEa 

A=D-V v F(D 1 ' v t oa ), (4) 

where v is the critical exponent of the correlation length 
for the 3D Ising universality class HI El El, and F(X) 
a crossover scaling function with X = D 1 l v t oa . Eq.QJ 
may qualitatively be interpreted using the finite size scal- 
ing principle ElEHsHs3 m which the film thickness D 
scales with the correlation length £ = £t^ l/ , where £ is an- 
other critical amplitude. To recover Eq.© from Eq.(0}, 
one notes that the scaling function F(X) must behave as 
F(X) oc X 13 for X — > oo. At a fixed small value of too, 
the D-dependence then cancels out from the equation, 
as it should. On the other hand, Eq.Q may also be re- 
covered from Eq.Q, by postulating that for small X a 
singularity occurs when X approaches X cr it, namely 

F(X) = f(X ~ X crit ) ft , X - A crit < 1, (5) 

with / another non-universal amplitude. This phe- 
nomenological assumption implies a scaling relation for 
the shift of the critical value of the polymer reservoir 
packing fraction 

A crit = D 1 '^ => = XcntD- 1 /". (6) 

Another scaling relation is implied for the critical am- 
plitude B(D), namely 

B(D) = fD^- V v . (7) 

It is clear that the crossover between both power laws, 
Eq.© and Eq.©, then also should occur when X is of 
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order unity, which implies very small values of al- 
ready when D is large. The region of — i<x, cnt where 
Eq.© and Eq.JHJ) then hold is extremely small. More- 
over, the general experience with problems of this kind 
is that a crossover never occurs abruptly HI |H |H H| , 
but rather spans several decades of the corresponding 
crossover scaling variable (X — X cr [ t in our case). In 
fact, if not a large enough range of this crossover scaling 
variable is accessible, one will instead observe a power 
law with "effective exponents" and "effective critical am- 
plitudes" , 

A(D) « 4ffi£ cff (8) 

where (3 2 < Pes < P pHj- The effective exponents do 
not have a fundamental deep meaning, of course, since 
their values depend on the range of tjj that is used for 
the analysis in terms of Eq. (JHJ) , and hence are not really 
defined unambiguously. 

III. MODEL AND SIMULATION METHOD 

We consider a mixture of hard-sphere colloids with 
diameter <j c and effective polymer spheres with diam- 
eter of gyration cr p in the grand canonical ensemble. 
Throughout this work, the colloid diameter a c is taken 
to be the unit of length. In the grand canonical ensem- 
ble, the volume V and the respective fugacities, z c and 
z p , of colloids and polymers are fixed, while the num- 
ber of particles in the system fluctuates. Following con- 
vention, the polymer fugacity is expressed in terms of 
a related quantity called the polymer reservoir packing 
fraction if p = 7rZper p /6. We also introduce the colloid 
packing fraction r/ c = ira^N c /(6V), with 7V C the num- 
ber of colloids in the system. The particles interact 
via potentials that were originally proposed by Asakura 
and Oosawa [65| (AO) and later, independently, also by 
Vrij (6(| . In this description, the so-called AO model, 
hard-sphere interactions are assumed between colloid- 
colloid and colloid-polymer pairs, while polymer-polymer 
pairs can interpenetrate freely. The interactions are thus 
strictly athermal such that the temperature plays no role. 
Instead, in the AO model, the analogue of (inverse) tem- 
perature is played by the polymer reservoir packing frac- 
tion rjp. As is well known, at the coexistence colloid 
fugacity and for sufficiently large colloid-to-polymer size 
ratios q = a p /a Cl the AO model exhibits a phase sepa- 
ration into a colloid-rich phase (the colloidal liquid) and 
colloid-poor phase (the colloidal vapor), provided rf p ex- 
ceeds a critical value. Grand canonical Monte Carlo sim- 
ulations are well suited to study this transition, and when 
combined with finite size scaling techniques these simu- 
lations also allow for investigations close to the critical 
point. Recently, this approach was applied to the bulk 
AO model, i.e. in the absence of walls, and the critical 
point was determined for q — 0.8, as well as the univer- 
sality class, which was shown to be that of the 3D Ising 
model [Hill. 



In this work, grand canonical Monte Carlo simulations 
are used to study the AO model in confinement. To this 
end, we use a simulation box of dimensions L x x L y x L z , 
with L x = L y = L and L z = D\ the system volume thus 
equals V = DL 2 . To capture the effect of confinement, 
we implement a so-called "sandwich" or "thin film" ge- 
ometry. Here, periodic boundary conditions are applied 
in the x and y directions, while in the remaining z di- 
rection we place two parallel walls: one in the z = 
plane, and in the z = L z plane. This closely resembles 
Ref. |4y where capillary condensation and evaporation of 
the AO model are investigated. Compared to the bulk 
AO model, one additional parameter is thus introduced, 
namely the film thickness D. For a film with thickness 
D, the thermodynamic limit is defined as the limit where 
the lateral dimensions L x = L y = L of the film are taken 
to infinity. Throughout this work, the walls are taken to 
be hard, i.e. colloid-wall and polymer-wall overlaps are 
strictly forbidden. Note that this implies a strong at- 
traction between the colloids and the walls due to the 
depletion effect 46] . The simulation method of Ref. |67] is 
then applied to the confined system; the main ingredients 
are a grand canonical cluster move [6jj and a reweighting 
scheme l85l. 



IV. RESULTS 
A. Binodal 

For q = 0.8 and D = 5, the grand canonical distribu- 
tion PL{j] c \rf v , z c ) is measured, defined as the probability 
of observing a system with colloid packing fraction ?y c , at 
"inverse temperature" rf p and colloid fugacity z c . There 
will generally be finite size effects in the lateral dimen- 
sions L x = L y = L of the simulation box, denoted by 
the subscript L. Phase coexistence is established by tun- 
ing z c such that Pl (Vc lWpi z c ) becomes bimodal with two 
peaks of equal area 86] . The respective packing fractions 
rjl and r?' , of the colloidal vapor and liquid phase, are ob- 
tained from the average peak positions. Typical distribu- 
tions are shown in Fig.^ plotted as kBT\nPL{r] c \r]p, z c ), 
with fee the Boltzmann constant and T the temperature. 
In this way, the distributions correspond to minus the 
free energy of the system. The height Fl of the peaks in 
k^ThiP^ri^rfp, z c ), measured with respect to the mini- 
mum in between the peaks (arrow in Fig.^), thus reflects 
the free energy barrier separating the colloidal vapor from 
the liquid phase |87| . In the two-phase region away from 
the critical point, the peaks in Pi(ry c |?7 p , z c ) are well sepa- 
rated and the barrier Fl will be large. Upon approach of 
the critical point, by lowering ry p , the peaks move closer 
together and the corresponding barrier Fl decreases pro- 
foundly. 

To obtain the binodal, Pl (tj c |^ p , z c ) is measured for a 
range of ry p and the average peak positions are recorded. 
The result is shown in Fig. [21 For comparison, the bin- 
odal of the bulk AO model is also shown, together with 
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FIG. 1: Coexistence distributions of the confined AO model 
with colloid-to-polymer size ratio q = 0.8, lateral dimension 
L = 15, and film thickness D = 5 for several values of rf p as 
indicated. Note that we have plotted the natural logarithm of 
PL{jlc\ifp, Zc), multiplied by k^T. In the above distributions, 
the colloid fugacity is tuned in order to obey the equal- weight 
prescription [86|. The barrier Fl in the distribution corre- 
sponding to rf v = 1.0 marks the average height of the peaks 
with respect to the minimum in between the peaks. 



the bulk critical point taken from previous work [671 168| . 
Note the familiar finite-size rounding in the simulation 
data close to the critical point (to describe the binodal 
correctly near the critical point requires finite size scal- 
ing, which we postpone to Fig. I10|l . The phase diagrams 
in Fig. |21 reveal the typical behavior of a fluid confined 
between two plates that undergoes a demixing transition: 
The critical point shows a significant shift which is ac- 
companied by an inward shift of the binodal with respect 
to the bulk 111. 



B. Cumulant analysis 

To obtain the critical polymer reservoir packing frac- 
tion rf p cr of the confined system, the fourth order cu- 
mulant — (to 2 ) 2 / (m 4 ) is measured, with m — 
Vc — (Vc), a s function of rj p for various lateral dimen- 
sions L. The cumulant is obtained by taking appro- 
priate moments of the distribution PL{r] c \rf p , z c ). For 
example, the average colloid packing fraction may be 
written as {r) c )(L,rf p , z c ) = J °° r] c P L (r] c \r] p , z c ) drj Cl and 
similarly for the p-th order moment (m p )(L,i] p , z c ) — 

Note that the outcome 
and 
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Jo [Vc ~ (Vc)} P P~L(Vc\Vp, zc) dr), 
will generally depend on r/ p , the colloid fugacity z, 
the lateral system size L. 



FIG. 2: Binodal of the AO model with q = 0.8 in bulk 
and confinement. Open circles show the bulk binodal, where 
the black triangle marks the corresponding critical point 
(ric = 0.134; rf p = 0.766) obtained using finite size scaling 
[671 In8l |. Open and closed squares show the binodal of the 
confined system with film thickness D = 5 for two lateral 
dimensions L. The horizontal line marks the critical poly- 
mer reservoir packing fraction rf v cr of the confined system 
in the limit L —* oo obtained using the cumulant intersection 
method, see Fig. [3] Lines connecting the points serve to guide 
the eye. Curve m shows the coexistence diameter (r/]. + rjZ)/2 
of the confined system with L = 20. 



Fig. [3] shows the cumulant as function of rf p for sev- 
eral system sizes L. We emphasize that the measure- 
ments were taken using the colloid fugacity at which the 
equal- weight prescription [8(1 was obeyed. At the criti- 
cal point, the cumulant becomes system-size independent 
[88|. The intersections in Fig. [3| thus provide an estimate 
of the critical polymer reservoir packing fraction. We ob- 
tain rf p cr = 0.9238 ± 0.0010, where the error reflects the 
scatter in the various intersection points. This estimate 
is also shown in Fig. (horizontal line) . Defining the co- 
existence diameter as (77' + rj^)/2 and ignoring finite size 
effects in this quantity for the moment, the intersection 
of the horizontal line with curve m (marked with a cross 
in Fig.[2J) yields an estimate of the critical colloid packing 
fraction ^ CiCr « 0.159. Compared to the bulk system, the 
critical colloid packing fraction has shifted to a slightly 
larger value. 

The cumulant plot also provides evidence for the 
crossover scenario discussed in the introduction. From 
Fig. [3 we obtain U4 « 0.795 at the critical point, which 
is between the 2D and 3D Ising values, see Tabled In the 
limit L — > 00, U4 is expected to approach the 2D Ising 
value. However, in the (still moderate) system sizes ac- 
cessible in our simulations, "effective" critical behavior 
is observed instead, with properties between those of the 
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FIG. 3: Cumulant analysis of the confined AO model with 
q = 0.8 and film thickness D = 5. Shown is the fourth order 
cumulant f7 4 as function of rjL for several lateral dimensions 
L. The intersection point yields an estimate of r?p jCr . The 
inset shows the cumulant slope Y\ at the intersection point as 
function of the lateral dimension L. 



2D and 3D Ising universality classes. Additional confir- 
mation of the crossover may be obtained from an analysis 
of the cumulant slope Y\ = dU^/drf evaluated at the crit- 
ical value of rip. It is expected that Y\ oc L 1 ^, with v 
the critical exponent of the correlation length and L the 
lateral system size. To this end we have plotted, in the 
inset of Fig. |3 the cumulant slope Y\ at the critical value 
rf p cr = 0.9238 obtained above, versus the system size L. 
By performing a fit to the data, the exponent is measured 
to be f c ff ~ 0.93, which is already surprisingly close to 
the 2D Ising value. Obviously, u c s must be interpreted 
as an "effective" critical exponent. 



TABLE I: Critical exponents of the order parameter (/3), 
correlation length (u), and specific heat (a) for the two- 
dimensional (2D) and three-dimensional (3D) Ising model, as 
well as the mean-field values (MF) . Also listed is the value of 
the fourth order cumulant (U4) at criticality for the 2D and 
3D Ising model. 
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FIG. 4: Finite size dependence of the free energy barrier Fl 
between the colloidal vapor and liquid phase, for the confined 
AO model with q — 0.8 and film thickness D = 5. Shown 
is Fl as function of 1/L at the indicated value of rf p , with L 
the lateral dimension of the simulation box. The barrier was 
extracted from equal-area [86fl distributions Pl{Vc\Vp> z c), see 
also Fig. Q 



C. Free energy barrier and interfacial tension 

Next, we consider the free energy barrier Fl between 
the colloidal vapor and liquid phase. At the critical point, 
the grand canonical distribution scales with the system 
size L as jHlil 

Pl{i lc ) = b Q L^ v V\b Q L^ v r) c ), (9) 

where P^(rj c ) is the distribution PL(rj c \rf v , z c ) measured 
in the finite system at the critical values of rf p and z c , bo 
some non- universal constant, and V° a function indepen- 
dent of system size (in the present case of confinement, 
the critical exponents (3 and v assume 2D Ising values). 
Recall from Fig. ^that Fl is given by the peak-to- valley 
height in the logarithm of PL{rj c \rf v , z c ). The scaling form 
of Eq. @ thus implies that Fl becomes system-size inde- 
pendent at the critical point, providing an alternative 
route to locate 77^ cr (see Ref. |93| where this approach 
is applied to the Lennard-Jones fluid). To locate cr , 
the barrier is recorded as function of L for several val- 
ues of ?/p. At the critical value of a plateau should 
be visible. For the confined AO model, the latter is ver- 
ified in Fig. QJ which shows Fl as function of 1/L for 
various rj^ around the critical region. The figure shows 
an increase in Fl with system size at high rf p , and a 
decrease at low rf v . The plateau occurs in the interval 
rf p = 0.923 — 0.927. Although not very precise, this esti- 
mate is consistent with the previous result based on the 
intersection of the cumulant. 
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FIG. 5: main frame: Interfacial tension 7^ of the confined 
AO model with q = 0.8 and D — 5 as function of rf p (solid 
curve), as well as the corresponding bulk interfacial tension 
(dashed curve), inset: 700 of the confined system as function 
of the relative distance from the critical point t, on double 
logarithmic scales, where r7p iCr = 0.9241 in t was used. The 
dashed lines illustrate 2D and 3D Ising exponents. 



A more precise estimate of rf v cr may be obtained from 
the critical behavior of the interfacial tension 700 in the 
thermodynamic limit. Upon approach of the critical 
point, starting in the two-phase region, the interfacial 
tension is expected to vanish as 



7oo = IV\ n 



(10) 



with t = rfp/rfp cr — 1 the relative distance from the crit- 
ical point, critical amplitude To, and critical exponents 
listed in Table [I] For 2D Ising systems we thus obtain 
/i2D = 1 and for 3D Ising systems /_i3D — 2v, where in the 
latter the hyperscaling relation 2 — a = dv was used (with 
d the spatial dimension). For 3D bulk, the expected ex- 
ponent 2v was already confirmed by us in Ref. 1681 In 
the present case of confinement, however, the crossover 
scaling scenario implies a transition in the critical behav- 
ior of 700 from singular (/i = 2v) in three dimensions to 
purely regular (/j = 1) in two dimensions. This particu- 
larly affects the slope of 700 versus t at the critical point, 
which should be zero in 3D, and finite positive in 2D. 

In order to test if evidence for this change in critical 
behavior is present in our simulation data, we use the 
free energy barrier Fl to measure the interfacial tension. 
Following Ref. Igj], the interfacial tension in a con- 
fined system of thickness D and finite lateral dimension 
L 3> D equals jl = Fl/(2LD), where the factor of two 
stems from the use of periodic boundary conditions. The 
thermodynamic limit L — > 00 can be evaluated through 
an elimination of finite size effects using the extrapolation 
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0,2 

LD ' 



(11) 



where the constants a± and et2 are a priori unknown. The 
interfacial tension in the thermodynamic limit 700 (f?p) 
at a given value of ryp may thus be obtained by fitting 
Eci. (|ll|l to corresponding measurements of 7l(t7 p ) in fi- 
nite systems. We have applied this approach using differ- 
ent system sizes between L = 12.5 — 22 and furthermore 
assuming 0,2 = in Eq. (f 1 1|) . The latter choice is based on 
empirical findings that the logarithmic term in Ea. Hlljl 
is typically rather weak, at least for Ising-like systems 
|68l l94| . The result is summarized in Fig. [S] The main 
frame shows the thermodynamic limit interfacial tension 
7oo as function of rf p for the confined AO model, as well 
as the bulk tension taken from previous work 68]. Note 
the pronounced decrease in the interfacial tension of the 
confined system with respect to the bulk, a direct con- 
sequence of the upward shift in 7y pcr . For the confined 
system, the vanishing of 700 at the critical point yields 
?7p cr = 0.9241 which is consistent with our previous es- 
timate based on the intersection of the cumulant. More- 
over, the interfacial tension seems to vanish with finite 
slope, a point further emphasized in the inset where 700 
is plotted as function of the relative distance from the 
critical point. Note that the interfacial tension in the 
confined system is already well described by the 2D Ising 
exponent. 



D. Order parameter 

Another consequence of the crossover scenario is that 
the binodal should appear "flatter" , since the critical ex- 
ponent (3 of the order parameter for the 2D Ising model 
is smaller than in 3D. In this section, we use the finite 
size scaling algorithm of Kim, Fisher, and Luijten j^| 
to study the critical behavior of the order parameter in 
confinement. The algorithm is based on the dependence 
of the cumulant U4 on the temperature- like variable rf pl 
the colloid fugacity z c , and the system size L (recall that 
U4 is defined in Section TlV B(l . For fixed 7y p and L, it 
is straightforward to measure U4 and (?7 C ) as function of 
z c . A plot of U4 versus (r/ c ), which is thus parameter- 
ized by z c , reveals two minima separated by a maximum, 
see Fig. El The location of the minimum at low colloid 
packing fraction is denoted 77^ (L, r/^), with Q~ (L, r/^) the 
corresponding minimum value. Similarly, the location of 
the minimum at high colloid packing fraction is denoted 
ri+(L,rf p ), with Q + (L,rf p ) the corresponding minimum 
value. Note that the location of the minima, and the 
corresponding minimum values, depend on rf v and L, 
but obviously not on z c . In the thermodynamic limit 
L — > 00, U4 approaches 1/3 in the one-phase region (hor- 
izontal dashed lines in Fig. EJ) . On the phase-boundary, 
7y ( T(L,?7 p ) and 77+ (i, 77^) approach the thermodynamic 
values ??c(?7p) and vlivl,)! respectively, while Q~(L,rf ) 
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FIG. 6: Cumulant ratio U$ as function of the average col- 
loid packing fraction (fj c )- The data were obtained for the 
confined AO model with q — 0.8 at -q p = 0.93, film thickness 
D = 5, and several lateral dimensions L as indicated. Dashed 
horizontal lines correspond to the limiting value U4, = 1/3, see 
details in text. The inset shows the region around r/~(L,ri p ) 
on an expanded scale. 
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FIG. 7: Order parameter A as function of rf p for the confined 
AO model with q = 0.8 and film thickness D — 5. The main 
frame shows A in the thermodynamic limit as function of rf v 
on linear scales, where the curve through the simulation data 
is a fit to Eq. JSJ . The inset shows the result as function of the 
relative distance from the critical point t = rf p /rf p cr — 1, on 
double logarithmic scales, where the slopes of the lines reflect 
2D and 3D Ising critical exponents. 



and Q + {L,Tf v ) approach zero [9g. The L-dependence 
shown in Figgis consistent with this scenario, see also 
the inset. 

In order for the scaling algorithm to succeed, simu- 
lation data for at least three different system sizes are 
required. In this work, L = 15, 17.5 and 20 are used. In 
addition, measurements over a rather broad range in 77I 
are required, starting in the two-phase region and stretch- 
ing toward the critical point. Here, five different rf v are 
simulated per system size, evenly distributed over the 
range ?7p « 0.9 — 1.0. Estimates of properties at inter- 
mediate ?7p are obtained using histogram extrapolation 
|97|. The purpose of the scaling algorithm is to evaluate 
the order parameter A as function of rf p in the thermo- 
dynamic limit 



A(<) = lim 



Following Ref. [95|, we define the quantities 

?+(L,rf) + Q-(i,rf) 



Qmin{L, Tip) 

x (L,r) T p ) = Qmin(£, ?7p) In 



2A(?7p) 



(12) 

(13) 
(14) 

(15) 



The algorithm starts in the two-phase region, with a 
value of ?7p significantly above its critical value. The 
peaks in P£,(?7 c |f7p, z c ) are then well separated and the free 
energy barrier Fl will be large (see for example the distri- 
bution corresponding to rj p = 1.0 in Fig.^J. This regime 
is called the "gaussian limit" because PL{j] c \rfp, z c ) may 
be described by a sum of two gaussians in this case |95| . 
In the gaussian limit, it can be shown rigorously that 
the points (x, y) of different system sizes L, should all 
collapse onto the line y = 1 + x/2. Recall that A(?7 p ) 
in Ea. l|15|l is the order parameter in the thermodynamic 
limit at the considered 77^, precisely the quantity of in- 
terest, which may thus be obtained by fitting until the 
best collapse onto y = 1 + x/2 occurs. Next, rfp is chosen 
closer to the critical point, the points (x, y) are calculated 
as before, but this time around A(?? p ) is chosen such that 
the new data set joins smoothly with the previous one, 
yielding an estimate of the order parameter at the new 
rL. This procedure is repeated as closely as possible to 
the critical point, where A vanishes, yielding an estimate 

of »7p,cr- 

For the confined AO model, the output of the scal- 
ing algorithm is illustrated in Fig. [7| The main frame 
shows the order parameter as function of 7f p on linear 
scales. The dashed curve is a fit to the simulation data 
using Eq.@ from which ifp cr = 0.9223, B cS = 0.173, and 
/3 c ff = 0.17 are obtained. As before, f} c g plays the role of 
an effective critical exponent. Note that (3 c g is already 
rather close to the pure 2D Ising value. This point is 
emphasized in the inset of Fig. which shows the order 
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parameter as function of the relative distance from the 
critical point t, where rj^ cr = 0.9223 in t was used. Also 
included are power laws illustrating 2D and 3D Ising crit- 
ical exponents. As expected, the simulation data slowly 
approach the slope of the 2D Ising exponent. By perform- 
ing additional simulations using larger lateral dimensions 
L, the data can in principle be extended to smaller val- 
ues of t, where the pure 2D Ising exponent will become 
visible. However, such simulations are computationally 
very demanding and beyond the scope of the present in- 
vestigation. In contrast, adding data at larger values of t 
in order to observe the 3D Ising exponent is not possible, 
since then we leave the critical regime. In hindsight, the 
thickness D = 5 considered here is too small to observe 
the full crossover from 3D to 2D Ising critical behavior. 
For such thin films, the critical behavior is essentially 2D 
Ising. The crossover scaling is expected to be visible only 
in much thicker films, where 2D Ising behavior shows up 
at extremely small t. 



In addition to the order parameter, the scaling algo- 
rithm also yields y as function of x. The latter func- 
tion, or scaling curve, is significant because it is uni- 
versal within a universality class. For bulk 3D fluids, 
belonging to the 3D Ising universality class, universality 
of the scaling curve has been verified for the hard-core 
square- well (HCSW) fluid [95j, the restricted primitive 
model (RPM) [95|, the decorated lattice gas [95j], the 
AO model [9j| , and the Widom-Rowlinson mixture [9^ • 
In the present case of confinement, however, the scaling 
curve is expected to deviate profoundly from the bulk 
3D Ising form. The latter is verified in Fig. |SJ which 
shows the scaling curve of the confined AO model ob- 
tained in this work, to geth er with the scaling curve of the 
3D bulk HCSW fluid [9^. Following the convention of 
Ref. I95I the scaling curve has been raised to an exponent 
—(f) — —1/(3, with (3 = 1/8 the critical exponent of the 
2D Ising model. An important feature of Fig. [H] is that, 
for small x, the data of the confined AO model correctly 
approach the limiting form y — 1 + x/2. In addition, 
we observe that the data from the three different sys- 
tem sizes have collapsed accurately onto a single curve. 
As expected, the scaling curve of the 3D bulk HCSW 
fluid differs profoundly from the one of the confined AO 
model, a direct consequence of the different universality 
classes. Note in particular the large difference in x c at 
which the scaling curve vanishes. For the HCSW fluid, 
x c 0.286 95], which exceeds the value of the confined 
AO model x c w 0.165 by over 60%. For the 2D Ising 
model, x c « 0.46 is reported [95|, which overestimates 
our value significantly. Note, however, that the scaling 
curve of Fig. [S] must be regarded as an "effective" scaling 
curve, and that deviations from the pure 2D Ising form 
are to be expected. 
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confined AO 
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3D bulk HCSW 
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FIG. 8: Order parameter scaling curve y - ^, with (f> = 1/(3 
and = 1/8, for both the 3D bulk HCSW fluid [H, and the 
confined AO model (q — 0.8 and D — 5) of this work. Also 
shown is the exact small x form y = 1 + x/2 of the gaussian 
limit. 



E. Coexistence diameter 

Finally, we turn to the critical behavior of the coexis- 
tence diameter S = (rj l c + rj^ )/2, in general given by |100| 

S = Vc, cr (1 + A 2p t 2 P + A^t 1 -* + A x t) , (16) 

with r/ c cr the colloid packing fraction at the critical point, 
t the relative distance from the critical point, and non- 
universal amplitudes A{. The term proportional to t 2 @ 
is due to pressure mixing, and for systems where pres- 
sure mixing is absent A 2 p — 0. It is not yet clear which 
features determine the degree of pressure mixing in a 
fluid. Of the bulk 3D fluids where this issue has been 
investigated, only the RPM exhibits substantial pressure 
mixing |9^. In the decorated lattice gas, pressure mix- 
ing is absent [9^|, and the same seems to be the case for 
the Widom-Rowlinson mixture j9^. Simulations of the 
HCSW fluid [13 and the AO model ^ point to rather 
weak pressure mixing. 

To obtain the coexistence diameter of a bulk 3D fluid 
is still challenging. In the present case of confinement the 
situation is even more subtle. Assuming negligible pres- 
sure mixing, the critical behavior of the diameter is dom- 
inated by t 1_Q . The crossover scaling scenario then im- 
plies a transition from weak singular behavior a « 0.109 
in 3D, to purely regular behavior a — in 2D, see TableQ] 
On the other hand, if pressure mixing is important, the 
diameter remains singular and dominated by t 2 ^ , with ul- 
timately P = P2 the critical order parameter exponent of 
the 2D Ising model. To determine which of these scenar- 
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FIG. 9: Coexistence diameter Sl given by Eq.(|THJ of the 
confined AO model with q = 0.8 and D = 5 obtained in finite 
systems of lateral dimension L as indicated. The arrow marks 
the estimate 77p. cr « 0.9223 obtained from Fig. |7| 



ios is realized in the confined AO model, we again use the 
finite size scaling approach of Kim, Fisher, and Luijten 
[95). Note that, in addition to the order parameter, these 
authors also present an algorithm to extract the diame- 
ter. As before, the algorithm generates a scaling curve, 
starting with data obtained well away from the critical 
point, and then recursively working its way down toward 
criticality. However, the quantities needed to construct 
the scaling curve are different. In particular, they involve 
the asymmetry factor 



A min (L,77p) 



Q+(L,Tf p )-Q-(L,Tf p ) 
Q+{L,rQ + Q-{L,%) 



(17) 



with Q^(L, r)p) defined previously. Unfortunately, for the 
confined AO model, we were unable to extract the diam- 
eter in this way. Closer inspection of our data revealed 
that A m i n as function of ?y p changes sign from negative 
to positive upon approach of the critical point, and this 
makes the procedure numerically unstable. In contrast, 
for 3D bulk systems, A m i n remains positive (at least for 
the AO model and the Widom-Rowlinson mixture) and 
so the problem does not occur there. 

Hence, our attempt to extract the critical behavior of 
the diameter in confinement is unsuccessfu l. Ins tead, we 
follow the more pragmatic approach of Ref . llOll and sim- 
ply show in Fig. El the diameter of the finite system 



(18) 



with rf^(L,rf v ) defined as before. Obviously, these data 
cannot be used to extract (effective) critical exponents, 
but the value r] c . cr w 0.159 of Section TlV Bl derived from 
equal-area distributions (77^77^, z c ) seems confirmed. 



_ Q. 




FIG. 10: Thermodynamic limit binodals, obtained using fi- 
nite size scaling, of the AO model with q — 0.8 in bulk and 
confinement. The dashed curve shows the bulk binodal, the 
solid curve is the binodal in the confined system with thick- 
ness D = 5. Triangles mark the corresponding critical points; 
squares and circles are raw simulation data obtained in finite 
systems away from the critical point. 



V. DISCUSSION AND CONCLUSION 

In this work, we have investigated the critical behav- 
ior of a colloid-polymer mixture confined to a thin film of 
thickness D = 5. The main finding is that for such thin 
films, 2D Ising universality is clearly visible. The latter 
is manifested by the critical exponents of the correlation 
length, the interfacial tension, and the order parameter. 
Since the order parameter exponent in 2D Ising systems 
is smaller than in 3D, the binodal in the confined sys- 
tem should appear "flatter". To emphasize this point, 
we have combined the order parameter data of Fig. [7| 
with the coexistence diameter of the largest system in 
Fig. and constructed the binodal in the thermody- 
namic limit. The result is shown in Fig. 1101 as the solid 
curve, where the closed triangle marks the location of the 
critical point. As expected, away from the critical point, 
the binodal obtained via finite size scaling joins smoothly 
with the raw finite-size simulation data (squares). For 
comparison, the binodal of the bulk system in the ther- 
modynamic limit is also shown, taken from previous work 
[68j . The different curvature of the binodals should be 
detectable in experiments. The result of Fig. may 
furthermore be relevant to Gibbs ensemble simulations 
of fluids in confined geometry. Here, the critical point is 
typically determined via a fitting procedure assuming 3D 
Ising exponents. In contrast, Fig. indicates that for 
thin films, extrapolations using 2D Ising exponents may 
be more appropriate. 
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In addition to a flatter binodal, the location of the 
critical point also changes with respect to the bulk. Our 
results indicate a pronounced shift of the critical point 
toward higher values of rf p , as well as a slight increase 
in the critical colloid packing fraction ?7 C!Cr . The increase 
in rf v is consistent with previous work on capillary con- 
densation in the AO model that was based on DFT and 
Gibbs ensemble simulations [46|. The behavior of ?7 ClCr 
is more subtle. For films with D > 5cx c , DFT shows an 
increase in ^ c cr with respect to the bulk value, while for 
very thin films a decrease is predicted |46j| . It is not obvi- 
ous if the corresponding Gibbs ensemble simulations also 
follow this trend [46|. At first sight, the increase in ?7 C!Cr 
observed in our simulations seems consistent with the 
trend predicted by DFT for films with D > 5er c . How- 
ever, for very thick films, f] c cr must ultimately approach 
its bulk value again, and it is not clear how the DFT of 
Ref. l4fil approaches this limit. Interestingly, recent Gibbs 
ensemble simulations of the confined Lennard-Jones fluid 
with film thickness D = 12 particle diame ters, show a 
pronounced decrease in the critical density 1 1 1| - which 
differs from both the present simulation result, and the 
DFT trend for this rather thick film. 

All in all, the critical density seems to depend quite 
sensitively on the details of the particle and wall interac- 
tions, as well as on the film thickness. Key to a reliable 
estimate of the critical density is a precise description of 
the coexistence diameter. The latter may be obtained 
using the finite size scaling approach of Ref. |9a as was 
recently demonstrated for 3D bulk fluids 95, 98, 99] . The 
issues pertaining to the shift of the critical density in con- 



finement inspired us, in Sectioi mV El to apply this scaling 
algorithm to the confined AO model. Unfortunately, we 
were unable to extract the diameter because the scaling 
algorithm of Ref. seems to behave profoundly different 
in confinement; the source is a numerical instability aris- 
ing from a change in sign of the asymmetry factor given 
by Ea. (|17(l . At this point, we see no reliable way to ex- 
tract the coexistence diameter of the confined AO model; 
nor do we understand the significance of the change in 
sign in A nlin . To resolve these issues would be the subject 
of further work. 

Regarding the order parameter, no difficulties were en- 
countered in applying Ref. |95| to the confined AO model. 
This was demonstrated by the accurate collapse of the 
data from different system sizes onto a single scaling 
curve. Also the estimated effective critical exponent p o g 
is consistent with the crossover scaling scenario. Never- 
theless, the large discrepancy in x c at which the scaling 
curve vanishes, between the confined AO model and the 
2D Ising model, is concerning. To understand the source 
of this discrepancy, too, would require further work. 
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